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Antideuterons are a potential messenger for dark matter annihilation or decay in our own galaxy, 
with very low backgrounds expected from astrophysical processes. The standard coalescence model 
of antideuteron formation, while simple to implement, is shown to be under considerable strain by 
recent data from the LHC. We suggest a new empirically based model, with only one free parameter, 
which is better able to cope with these data, and we explore the consequences of the model for dark 
matter searches. 
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I. INTRODUCTION 


The use of antideuterons for indirect detection of dark 
matter (DM) was first suggested in [T]. Despite the low 
yield per annihilation or decay, antideuterons can be an 
important discovery channel due to the extremely low 
astrophysical background. Presently, the AMS-02 exper¬ 
iment is taking data that can improve on current upper 
bounds for the antideuteron flux at the Earth. 

There are several uncertainties at play when calculat¬ 
ing the resulting bounds on dark matter models. The 
most significant is the uncertainty in propagation mod¬ 
els, while the second is the antideuteron formation model. 
The dark matter halo uncertainty can also be large. We 
will here concern ourselves with the formation model. 

The formation of antideuterons is commonly described 
using the so-called coalescence model, which harks back 
to the 1960s [H [3]. In this simple phenomenological 
model, any antiproton-antineutron pair with momen¬ 
tum difference \jfp — Pn\ < Po, will combine to form 
an antideuteron. The coalescence momentum pq, typ¬ 
ically evaluated in the centre-of-mass (COM) frame of 
the antinucleons, is a free parameter that must be hxed 
by calibration against experimental data. While modi¬ 
fied slightly over the years, the coalescence model is still 
state-of-the-art. In calibrating pq on (relatively) modern 
experimental data, it has been found that while all avail¬ 
able datasets can individually be consistently described 
by some value of Pq, there is no consistent pg between 
datasets Hig. This has been explained by differences 
in the event generators used to simulate the data, and 
by the different physical properties of the processes mea¬ 
sured, e.g. production from a colorless e“'"e“ initial state 
versus production in pp-scattering. 

In this work, we will show that new data on deuteron 
and antideuteron production from the ALICE experi¬ 
ment at the Large Hadron Collider (LHC) [7], cannot be 
well described by the coalescence model. Thus, for the 
first time, challenging the model in a single experiment. 
We will present a new, empirically based model that de¬ 
scribes (anti)deuteron formation as a probabilistic pro¬ 
cess, and show that it is capable of successfully describing 
the new ALICE data. This model will then be applied to 
make predictions on the antideuteron flux from a generic 


annihilating dark matter model, and we will compare it 
to the predictions of the coalescence model. Along the 
way we will also comment on the potential usefulness of 
future data on deuteron production in order to explore 
our model further. 

In Section|ll]we will begin by reviewing the coalescence 
model and some of its recent modifications. We then go 
on to describe the basis of our new model in Section Hill 
In Section [TVl we proceed by comparing the calibration 
of the two models on a selection of the available datasets. 
Section|V]describes the resulting cosmic ray antideuteron 
flux from dark matter in our model, comparing it to the 


coalescence model, before we conclude in Section VI 


II. THE COALESCENCE MODEL 

In its initial form, as it was first applied to deuteron 
production in heavy ion collisions, an additional assump¬ 
tion of isotropic and uncorrelated antiproton and an¬ 
tineutron spectra was used in the coalescence model to 
obtain an analytical expression for the antideuteron spec¬ 
trum in terms of the antiproton and antineutron spectra. 
These assumptions have, however, been show not to hold 
in processes relevant to indirect DM detection [5], and 
the coalescence condition should therefore be applied to 
ph-pairs on a per-event basis. 

As has been show in Refs 015 ], tuning Po against ex¬ 
periments measuring different collision processes at dif¬ 
fering energy scales does not give a consistent best fit 
value. Moreover, the coalescence model is sensitive to 
two-particle correlations for (anti)baryons, arising from 
the structure of the hadronization models in the Monte 
Carlo event generators used [S]- Tuning different event 
generators to the same experimental data will therefore 
typically not give the same best fit values for the coales¬ 
cence momentum. As discussed in Ref. |3], hadroniza¬ 
tion parameters in Monte Carlos are usually not tuned 
to measured two-particle correlations, where they exist, 
nor indeed with any specific emphasis on reproducing 
(anti)nucleon spectra. Tuning hadronization parameters 
specifically towards antideuteron production is therefore 
a prospective way of achieving better consistency in fits 
to experimental data, as well as better agreement be¬ 
tween different Monte Carlos. 
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It was pointed out by the authors of Ref. [S] , that spa¬ 
tial separation should also be taken into account when 
evaluating the coalescence condition. Nuclear interac¬ 
tions take place on scales of a few femtometers, while 
weakly decaying particles will typically have macroscopic 
decay lengths. Their decay products will therefore be 
produced too far from the primary vertex to have a 
chance of interacting with particles produced at the pri¬ 
mary vertex. For this reason, weakly decaying particles 
should be considered stable in the context of coalescence. 
As an alternative, the authors of Ref. [5] implement an 
explicit condition on the spatial separation between the 
antinucleons of Ar < 2 fm in their coalescence model, 
which in principle is a more correct approach. How¬ 
ever, since most Monte Carlos do not model the space- 
time structure resulting from showering and hadroniza- 
ton, and we expect very few antideuterons to be produced 
by decaying final states, we expect the two approaches to 
be more or less equivalent. 

The coalescence model for antideuteron production de¬ 
scribes a 2 —>■ 1 process, which does not preserve energy- 
momentum. This issue is not much discussed in the liter¬ 
ature, but is usually solved by requiring momentum con¬ 
servation, Pd = Pp+Pn, and calculating the antideuteron 

energy through = \J\Pd\‘^ + mg, implicitly assuming 
that the excess energy is somehow disposed of at a later 
point. A more satisfactory description is to consider this 
as a radiative capture process pfi —>■ dy, which is the 
dominating antideuteron formation process at the low 
COM momentum differences required by the coalescence 
model. For a full kinematical description, the magnitude 
and direction of the photon recoil must be taken into ac¬ 
count through four-momentum conservation. However, 
for antideuteron kinetic energies well above po, the effect 
is negligible. Any spin correlations in the COM system 
will affect the angular distributions of the final state par¬ 
ticles with respect to the initial state, and should in prin¬ 
ciple also be taken into account. However, we see no a 
priori reason for such a correlation, and the effect will be 
washed out in the lab-frame by the generally large boost. 


III. AN EMPIRICAL, CROSS SECTION BASED 
MODEL 

A. The model 

In the coalescence model, antideuteron formation is 
classically deterministic, and the probability that a pn- 
pair will form an antideuteron can be expressed as a step 
function in the COM momentum difference between the 
antineutron and antiproton, 

P{pn ^ d\k) = 9{po - k), (1) 

where k = \pp — Pn|cOM- From quantum mechanics, one 
would not expect a relation like this, but rather a for¬ 
mation probability that depends on the wave function 


overlap of the initial state nucleons, and varies as a func¬ 
tion of k, just as in an ordinary scattering process. We 
expect this probability to be proportional to the cross 
section for the corresponding capture process pfi dX, 

P{pn ^ dX \ k)(x apf,^dx{k). (2) 

As an alternative to the coalescence model, we there¬ 
fore propose a model in which the combination of a 
ph-pair with COM momentum difference k into an an¬ 
tideuteron is a random event with a probability given 
by 


P{pn ^dX \k) = ^ 


where (Tpn^dxi^) i® the sum of cross sections for pn- 
processes with an antideuteron in the final state, and 
the constant of proportionality (Tq i® a free parameter to 
be fixed through calibration against experimental data, 
analogous to po in the coalescence model.^ 

For low values of fc, the relevant process is the ra¬ 
diative capture process ph —>■ dy. For COM energies 
above the pion production threshold, instead processes 
with hadronic final states pfi —>■ d{NTT)'^ dominate. At 
these energies, antideuterons are actually more efficiently 
produced through pp and hh processes with d{NTT) final 
states, and these processes must therefore also be taken 
into account. The cross sections decrease with increasing 
number of final states, and experimental data also be¬ 
come significantly more sparse. In this work, we will as 
a result only consider the antideuteron production pro¬ 
cesses listed in Table HI 


1) pn —>■ dy 

2) pn —> Jtt® 

3) pn —> d'K^'K~ 

4) pn —>■ d7r°7r° 


5) pp —>■ d'K~ 

6 ) pp —> d'K~'K° 

7) fin —>■ diT^ 

8 ) fin —>■ d-K^n^ 


TABLE I: Processes considered in this work. 


For a given antinucleon pair, the probability that it will 
form an antideuteron though a process i from Table |l] is 
in our model given by 

P{NiN 2 dX, I k) = 

<Po 

where iVi and N 2 are the species of the two antinucle¬ 
ons, and Xi represents the other final state particles in 
the given process. The free normalization factor tJo is 
assumed to be the same for all processes. The energy of 
the produced antideuteron depends on the kinematics of 


^ While (TO should in principle be calculable, it will in practice 
depend on properties of the wave-functions of the incoming nu¬ 
cleons. 
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the relevant process, and we will discuss this separately 
for the different processes in the following sections. 

Little or no data is available on the antinucleon pro¬ 
cesses we consider here, and we will therefore be basing 
our model on fits to data on the charge conjugate pro¬ 
cesses under the assumption = '^NiN^^dx- 


B. The pn dy process 

We have found only a small amount of data on the 
pn —>■ dj process, and then only at low energies. This 
alone is not sufhcient to make a fit of the cross section 
as a function of k. However, for the inverse process of 
photodisintegration, dy —>■ pn, a large amount of data 
is available, and can be used through application of the 
principle of detailed balance — see e.g. Ref. [10] for a de¬ 
tailed description. The principle implies that given time 
reversal invariance of the interaction, the cross section 
for a process <7{Aa — )■ Bb) is related to the cross section 
for the inverse process through 

2 

a{Aa Bb) = Aa), (5) 

9Ada Pi 

where pi is the momentum, and pi is the number of spin 
states of particle i; for massive particles, pi = {2si -\- 1). 
All quantities are given in the COM frame. Cross sec¬ 
tions are invariant under Lorentz boosts along the beam 
axis, and most experimental cross sections can therefore 
be used at face value here. In the derivation of the above 
expression, applicability of perturbation theory is typ¬ 
ically assumed, but the principle can be shown to be 
valid also when perturbation theory breaks down, pro¬ 
vided that averages over all spin variables have been per¬ 
formed [TT| . 

Applying the principle to the process pn -A d'y, we 
have Sp = s„ = i and = 1, giving = g^ = 2 and 
Pd = 3. While the photon has spin 1, it is massless and 
thus only contributes two polarization states, pj = 2. 
Detailed balance then gives the relation 

3 p‘^ 

a{pn^ d'y) =-^aid'j ^ pn), ( 6 ) 

^ Pn 

which is frequently used in the literature on radiative 
capture and deuteron photodisintegration. 

Large amounts of experimental data on deuteron pho¬ 
todisintegration can be found in the literature, however, 
some of the experiments are in tension with each other. 
In order to be able to make a fit, it is necessary to prune 
the data down to a consistent dataset. Lacking better 
information, our approach is therefore to only use data 
from the most recent experiment in energy ranges where 
the experiments are in tension. For consistency, we dis¬ 
card the entire datasets from removed experiments, not 
only the points that are in tension with other experi¬ 
ments. Our final set of experimental data consists of 


Parameter 

Value 


a_i 

2.30346 


0.0 

-9.366346 

xlO^ 

Oi 

2.565390 

xlO^ 

a2 

-2.5594101 

xl0‘‘ 

03 

1.43513109 

xlO® 

a4 

-5.0357289 

xlO® 

05 

1.14924802 

xlO® 

QiQ 

-1.72368391 

xlO® 

Or 

1.67934876 

xlO® 

as 

-1.01988855 

xlO® 

Og 

3.4984035 

xlO® 

Oio 

-5.1662760 

xlO'^ 

bi 

-5.1885 


b2 

2.9196 



TABLE II: Best fit values to the parameters given in Eq. Q . 

radiative capture data from Refs. [T^HTH] and photodis¬ 
integration data from Refs. 

After applying the principle of detailed balance to 
translate the photodisintegration data into radiative cap¬ 
ture cross sections, we perform a least squares fit to the 
combined radiative capture data using the function 

<^np^d^{K) ^ I ZZ-1 : K < 1.28 . , 

(I/ib) \ exp(—biK — b 2 K^) :k>1.28, 

where n = k/{l GeV). We chose to use an exponential 
form above /c = 1.28 to ensure that the function does not 
unphysically diverge or obtain negative values at high 
energies. Due to the k~^ term, the fit function for the 
cross section clearly goes to infinity as k approaches 0. 
We therefore take care to restrict the antideuteron pro¬ 
duction probability to P{np -A d'y \ k) < 1 when using 
this fit in Eq. The best fit parameter values for our 
dataset can be found in Tab. |IT| and give an excellent fit 
of = 51.8 for 83 degrees of freedom. Note that since 
the fit was made to data spanning 6 orders of magnitude 
in energy, the parameter values are rather finely tuned, 
and must therefore be used at the given level of precision. 
The experimental data, as well as our fit are plotted as 
a function of k in Fig. [l] The peak in the cross section 
near I GeV is due to the delta-resonance — processes in 
which one of the nucleons is excited to a delta resonance, 
via virtual pion exchange, as seen in Fig. [^ 

In this model, in contrast to the coalescence model, an- 
tideuterons can be produced at values of k well into the 
GeV range, which leaves a potentially large amount of ex¬ 
cess energy to be radiated off by the photon. This in turn 
gives a sizeable recoil that must be taken into account by 


^ We were unable to reliably extract the errors from Ref. m, 
and instead assumed 5% errors, which are typical for similar 
experiments. 
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FIG. 1: Fit to experimental data for the deuteron radiative 
capture cross section as function of nucleon momentum dif¬ 
ference, k, in the COM frame. Circles show the photodisin¬ 
tegration data, while diamonds show radiative capture data. 



FIG. 2: Feynman diagram for the delta-resonance in radiative 
capture. Ni are here antinucleons. 


requiring four-momentum conservation. In application of 
the model we let the antideuteron and photon be emitted 
back-to-back in the COM system, in a random direction 
drawn from an isotropic distribution. 


C. N 1 N 2 —> dn processes 

The pn —)■ , nn —)■ dTT~ and pp —>■ d7r+ processes 

are related by isospin invariance through 


and 


pion masses. Very little data exists for the pn —>■ dn^ 
process, and we have not been able to find any data on 
the nn —>■ d'K~ process. A substantial amount of data 
is, however, available on the pp —>■ d7r+ reaction, and 
we will therefore use these data in combination with the 
above isospin relations to approximate the pn —>■ d7r° and 
nn —>■ d'K~ cross sections. The authors of Ref. [55] have 
already made a fit to the available data on the pp —>■ dTr"*' 
process, and we will here adopt their fit. They find the 
data to be well described by the function 


a{p) 


aif 

(c — exp(d? 7))2 -f e 


( 10 ) 


with the parameters given in Tab. Ill where p = g/TO.n.+ , 
and q is the momentum of the pion in the COM frame.^ 


Parameter Value 
a [pb] 170 

b 1.34 

c 1.77 

d 0.38 

e 0.096 


TABLE III: Best fit values from Ref. |26| to the parameters 
given in Eq. (101. 


The fit was made in the context of comparison to 
pn —>■ d7r° data, and was corrected for Coulomb repul¬ 
sion and phase space differences due to the differing pion 
and nucleon masses. These effects should in principle 
be re-applied to the pp-process, and an analogous phase 
space correction should also be applied when using the 
fit with the nn-process. However, these effects are only 
important near threshold for the process, and will effec¬ 
tively shift the threshold slightly in k. At high COM 
energies, the cross section is unchanged, and at the peak 
the corrections are only at the percent level. There is 
no reason to expect the (anti)deuteron spectrum to be 
sensitive to the precise position of the threshold, so for 
simplicity we will neglect these corrections here. We set 
the cross sections to zero below the kinematic thresholds 
for the processes, as this is not ensured by the fit. 

We have plotted the cross section fits for the processes 
as function of the COM momentum difference k in Figs.j^ 
and[^ The cross sections for these processes also peak 
at the delta resonance near k = 1 GeV, and, just as for 
the photon in the pn —>■ d'j case, the pion recoil must 
be taken into account. We again emit the antideuteron 
and pion back-to-back in a random, isotropically drawn 
direction in the COM frame, and determine the four- 
momenta from the kinematics. 


^nn^d'K~ ^pp^d'TT'^ ? 


( 9 ) 


see e.g. Ref. These relations are not exact, as the 
isospin symmetry is broken by the differing nucleon and 


^ When using the isospin relations, one should for consistency use 
in calculating rj also in the pn and nn d7r~ 

processes ESj. 
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D. A'^iA^2 —>■ dn-K processes 

For the N 1 N 2 —)■ dirTt processes, data are available 
on all but the nn —t d7r~Tr^ process. We here use 
pp d7r+7r° data from Refs. [221 [27H29], pn —)• dn~^'K~ 
data from Refs. [221 [211120]) and pn —>• dn^ir^ data from 
Refs. [201 Elj-^ There is unfortunately very little data 
available for ^/s > 2.5 GeV for all the processes. This 
makes fits to the pp —>■ dir^ir^ and pn —>■ dTr+Tr” pro¬ 
cesses particularly problematic. The exact locations and 
heights of the resonance peaks near ^/s = 2.5 GeV in 
these two processes are unclear, and for the pp —> d7r+7r° 
process, the lack of data at high energies makes the naive 
fit quite unstable. To improve on this, we again make use 
of isospin invariance. 

Isospin invariance predicts the relations [25j 

>-rf7r+7r“ f T '^^pp^d'K^•: (H) 


and 


xiTT TT^ >-d7r+7r° ? 


( 12 ) 


between the cross sections. Measurements of the pro¬ 
cesses in Eq. (Ill within the same experiment [22] have 


shown these cross sections to be quite sensitive to isospin 
breaking effects, leading to a ~ 25% deviation in this re¬ 
lation. If the isospin symmetry was exact, one could have 
used Eq. (Ill to make simultaneous fits to all three pro¬ 


cesses, but due to the isospin breaking we have not been 
able to obtain good fits in this manner. We therefore in¬ 
stead perform individual fits to each process, where we in¬ 


clude the data from the other processes through Eq. (11) 


for stability, but weighted down by a factor 1/100 in the 
X^. We have chosen the value of the weight to be large 
enough to guide the fits, giving reasonable high energy 
behaviour in the pp —> d7r"''7r*^ and pn —> dTr+Tr” chan¬ 
nels, but low enough to give good individual fits for the 
different processes. To further guide the fits, we also in¬ 
sert dummy data points at the kinematic cutoffs for the 
processes, with zero cross section, and errors of 1 /rb. 

We use the following functional forms for the fits, in¬ 
spired by [21], 


(j{k) = 


aK 


(c — exp(dA^))^ + e ’ 
for the pp —>■ d 7 T~^ 7 r^ and pn —>■ dir^ir^ processes, and 

^b 2 


(13) 


r{K) = 




02 K 


(ci — exp(diK))'^ + Cl (c 2 — exp(d 2 K))^ + 62 ’ 

(14) 


^ Many of the datasets are only available as plots, and for cases 
where errorbars cannot be resolved, we use the point size of the 
plot as an estimate for the error. 


for pn —d7r“^7r“, where again n = k/{l GeV). The 
best fit parameters for the different processes are listed 
in Tables [IVj |V| and |V^ The data points used and our 
fits to these points are plotted as functions of k in Fig. [^ 
No data are available on the nn —>■ dTT~'n^ process, and 
we are forced to make use of Eq. (121 to approximate 
the cross section for this process. As in the N 1 N 2 dir 
case, we set the cross sections to zero below the kinematic 
thresholds. 


Parameter Value 
a [pb] 2.855x10® 

b 1.311x10^ 

c 2.961x10® 

d 5.572x10“ 

e 1.461x10® 

TABLE IV: Best fit parameters for the pn —>■ d7r“7r“ process. 


Parameter Value 


ai [pb] 

6.465x10® 

bi 

1.051x10® 

Cl 

1.979x10® 

di 

5.363x10“ 

Cl 

6.045x10® 

02 [pb] 

2.549x10®® 

b 2 

1.657x10® 

C2 

2.330x10'^ 

d 2 

1.119x10® 

62 

2.868x10®® 


TABLE V: Best fit parameters for the pn —^ d7r'''7r process. 


Parameter Value 

a[ph] 5.099x10®® 

b 1.656x10® 

c 2.333x10'^ 

d 1.133x10® 

e 2.868x10®® 


TABLE VI: Best fit parameters for the pp —>■ dTr+Tr® process. 

As these processes have three-body final states, the 
kinematics become considerably more involved than in 
the previous cases. For a detailed review, we refer to the 
section on three-body decays in Ref. [22]- In processes 
with three-body final states, there can be angular correla¬ 
tions between the final states that depend on the matrix 
element for the process, and this is the case for the pro¬ 
cesses considered here. Dalitz plots from measurements 
of the pn and pp processes can be found for a few dif¬ 
ferent COM energies in Ref. [22], but these data are not 
sufficient to parameterize the deuteron GOM momentum 
distribution as function of energy. We therefore make 
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FIG. 3: Fits to cross section data. Left: pn —>■ d-K^Tv^, middle: pn —>■ d7r"'"7r , right: pp —> ci7r+7r°. 


the approximation of no angular correlations between the 
outgoing deuteron and pions, and draw the deuteron mo¬ 
mentum based on phase space alone. We determine the 
deuteron momentum by first drawing random invariant 
masses and uniformly within the kinematically 
allowed region. The momentum of the deuteron in the 
COM frame is then given by 

and we draw its direction from an isotropic distribution 
in the COM frame. 



FIG. 4: Fits to cross sections for np —>■ dX processes. 


E. Process contributions 

In the coalescence model, all antideuterons are by con¬ 
struction produced by ph-pairs with low COM momen¬ 
tum differences. Our model, on the other hand, has the 
majority of antideuterons produced close to the delta res¬ 
onance near k = 1 GeV. While radiative capture pn —^ d'y 
has a very high cross section at low values of k, the num¬ 
ber of available pn pairs drops very quickly for decreasing 



FIG. 5: Fits to cross sections for pp —>■ dX and nn —^ dX 
processes. 


values of k in the processes we have studied. This can 
be seen in Fig. where we show the number of possible 
antinucleon-antinucleon combinations in LEP events at 
the Z-peak, generated using Herwig++ 2.6.0, as func¬ 
tion of k and the combined momenta of the antinucleon 
pairs in the lab frame. The distribution peaks for values 
of k in the low GeV range and drops quickly for decreas¬ 
ing values of k. The result is that radiative capture at 
low values of k gives a very small contribution to the total 
antideuteron spectrum in the cross section based model. 
Instead, the peak in the number of antinucleon pairs is 
close to the delta resonance for all values of the combined 
lab frame momentum. Antideuteron production is thus 
dominated by the iViiV ^2 —>■ d{NTr) processes for more or 
less any antideuteron lab-frame momentum®. This holds 
true in all the experiments we consider in this work. 

Another notable feature in Fig. is that the total 
number of available ph-pairs is roughly a factor two larger 
than the corresponding numbers of pp and fin pairs, as 


® The total momentum of an antinucleon pair is an approximation 
for the momentum of the resulting antideuteron. 
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can be expected from pure combinatorics. In LEP events, 
antiprotons and antineutrons are produced with approx¬ 
imately equal probabilities. Picking two random antinu¬ 
cleons that each have equal probability of being an an¬ 
tiproton or antineutron is twice as likely to give a ph-pair 
than it is to give either a pp-pair or a fin-pair. Since the 
cross sections for antideuteron production in pp and fin 
processes are a factor two larger than the pn cross sec¬ 
tion in the delta resonance region, this implies that pn, 
pp and fin processes give similar contributions to the an¬ 
tideuteron spectrum in our model. 


F. Monte Carlo implementation 

The following is a step-by-step description of how our 
model should be applied to Monte Carlo events: 

• For each event, iterate over all possible unique 
antinucleon-antinucleon pairs, avoiding double 
counting pfp and hh-pairs. 

• For each antinucleon-pair in the event, calculate the 
momentum difference k between the antinucleons 
in their COM frame. Calculate the probabilities 
P{NiN 2 —f dXi I k) that the antinucleon pair will 
form an antideuteron for each relevant processes i 
listed in Table using Eq. 0. NiN'i^dXii^) is 
given by the fits in the previous sections,® and CTq 
has to be determined by fits against experimental 
data.^ 

• Draw a random number uniformly on the unit 
interval for each possible formation process. If 
Ti < P{NiN 2 —f dXi I k) for one of the processes, 
the antinucleon pair forms an antideuteron through 
that process. If < P{NiN 2 —f dXi \ k) for 
more than one process, pick one of the processes 
randomly using probabilities equal to the relative 
cross sections.® If an antideuteron is formed, ex¬ 
clude the involved antinucleons from being used in 
the formation of other antideuterons.® 

• Emit the antideuteron in a random, isotropically 
drawn direction in the COM frame. For two-body 


° Note that the cross sections for the processes with a single pion 
in the final state are parameterized on 77 = (where q is 

the pion momentum in the COM frame), rather than k. 

^ We will return to the fitted (jq values below. 

^ As the probability of having multiple successful processes is very 
low, the way in which a process is chosen in these cases has no 
significant effect on the final spectrum. 

^ The authors of m estimate that even at large values of po = 250 
MeV in the coalescence model multiple successful antideuteron 
candidates are found in less than 0 . 1 % of the events, and are thus 
a negligible problem. This also holds true in our model due to 
the low probabilities for each pair. It is interesting to speculate if 
these very rare events could be used to estimate the production 
of even heavier antinuclei such as ^He. 


final states, its energy and momentum is deter¬ 
mined by four-momentum conservation. For three- 
body final states, draw the antideuteron momen¬ 
tum randomly based on the available phase space, 
as discussed in Sec. EIdI 

• Boost the antideuteron to the lab frame. 


G. Extracting more information 


In the cross section based approach, the determination 
of whether or not an antinucleon-antinucleon pair will 
form an antideuteron is probabilistic. In a single event, 
there can be many possible iViV-combinations, each with 
a non-zero probability to produce an antideuteron. These 
probabilities will typically be very low, and in most cases, 
none of the combinations will produce antideuterons — 
the event has essentially gone to waste in the Monte Carlo 
statistics. Even in events where antideuterons are pro¬ 
duced, other W-combinations could also have been pos¬ 
sible, and this information would remain unused. This 
situation is similar to the one in the coalescence model, 
where extremely large statistics are needed in order to 
get a precise antideuteron spectrum. 

Moreover, for a given antinucleon pair, the energy of 
the resulting antideuteron is not fixed: since the COM 
frame of the antinucleon pair typically will be boosted, 
the antideuteron energy will be determined by the ran¬ 
domly chosen direction in which it is emitted. There 
is, in other words, much more information in each event 
than will be extracted by a single application of the an¬ 
tideuteron production model. 

In order to extract more information from these events, 
one can use weighted antideuteron events. For each 
Monte Carlo event we: 


• Set up a temporary histogram with the same bin¬ 
ning as the one used for the total antideuteron spec¬ 
trum (main histogram). 


• If the event contains more than one antinucleon, 
evaluate the event iVsamp times, following the pro¬ 
cedure in Sec. IIIF[ and adding any antideuterons 
to the temporary histogram. 


For a given bin, b, in the main histogram, the number of 
antideuterons can then be calculated as 


A'mc 

< = ^ (16) 
i=l 


where IVmc is the total number of Monte Carlo events. 


Wb,i = 


X ^a.mr 


(17) 


is the contribution to this bin from Monte Carlo event i, 
and is here the number of antideuterons in bin b in 
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FIG. 6: Number of available antinucleon pairs in e''"e“-collisions at the Z-resonance in Herwig++, as a fnnction of antinucleon 
momentum difference k in the COM frame and total momentum of the pairs in the lab frame. Left: ph-pairs, middle: p^pairs, 
right: hh-pairs. The plots have a shared normalization, and are normalized to a maximnm valne of 1. 


the temporary histogram of event i after iVsamp evalua¬ 
tions. The error on the number of antideuterons in bin b 
in the main histogram is then given by 




. ^MC 
\ i=i 


(18) 


This error is found under the assumption that IVsamp is 
large enough to give a representative sample of the an- 
tideuteron spectrum in each event, but in practice we 
find it to give a good error estimate even with a rela¬ 
tively low number of A^samp = 10. Using this method, we 
have found that it is possible to achieve the same level of 
precision with an order of magnitude fewer Monte Carlo 
events. This would not be possible in the coalescence 
model, as the antideuteron formation in that model is 
deterministic, and no more information could be gained 
by evaluating the same event multiple times. 


in their own experiment, and thus rejected the arguments 
for the introduction of the coalescence model. This con¬ 
troversy has apparently been more or less forgotten, and 
the coalescence model has remained state-of-the-art up 
to today. 

Experimental cross sections have also been used to 
estimate deuteron production in later works, such as 
Ref. [3D]. Here, the authors model the deuteron forma¬ 
tion in the np —> d7r''"7r“ process as a np —>■ NNtt process 
followed by a NN dir process. They use a conventional 
scattering model to describe the first step of the process, 
and then use experimentally measured cross sections for 
the NN —7^ dn processes to model the deuteron forma¬ 
tion. Using this model, they obtain results in reasonable 
agreement with the experimental measurements. 


IV. CALIBRATION AGAINST 
EXPERIMENTAL DATA 


H. An historical aside 

Modeling deuteron production based on experimen¬ 
tally measured nucleon-nucleon cross sections was dis¬ 
cussed in the original coalescence paper by Schwarzschild 
and Zupancic [3| from 1963, but they found this approach 
to yield too few antideuterons. They thus argued for 
the presence of a mechanism in which interactions with 
the surrounding nuclear matter affects the production 
of deuterons, and introduced the coalescence model as a 
phenomenological, simplified version of a model proposed 
by Butler and Pearson [33] . 

This approach was criticized by Kamal et al. in an 
article from 1966 |3S|. Here, they point out that the 
NN —>■ diT processes have a resonant behaviour, and 
that Schwarzschild and Zupancic had significantly un¬ 
derestimated the cross section for these processes. They 
also note the necessity of including contributions from 
pp —>■ dTr~^ and nn —>■ d'K~ processes. Taking the resonant 
behaviour and the extra processes into account, they ob¬ 
tained results in agreement with the deuteron production 


The coalescence model and the cross section based 
model each have a free parameter that needs to be tuned 
against experimental data. We here present the best fit 
parameter values for various experiments and two dif¬ 
ferent Monte Carlo event generators, giving necessary 
details on the experiments and event generation for re¬ 
producibility. We use the Herwig++ 2.6.0 |5D1I37] and 
Pythia 8.186 [38l|39] event generators with default set¬ 
tings, unless stated otherwise. In order to be able to 
compare the fits for the two models, we only use the ex¬ 
perimental uncertainty in calculating for the fits. This 
is to avoid bias in the due to differing statistics from 


After completion of this work, we have also been made aware of 
a paper by Gugelot and Paul Bol from 1993 that suggested a 
cross section based deuteron formation model similar to ours for 
use in Monte Carlos. At the time, there was unfortunately not 
sufficient experimental data available to properly test the model, 
and the work has largely gone unnoticed. We would like to thank 
Sebastian Wild for notifying us of this work. 
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the event generation. The statistical uncertainty is in all 
cases small relative to the experimental uncertainty, so 
the effect of this should be small. The resulting values 
for po and cto are listed in Tables |VII| and |VIII[ As the 
antideuteron spectrum in the cross section based model 
scales as oc 1/cto, we present the results for this model in 
terms of I/cto: rather than ag itself. 

As discussed in Sec. |Tlj antinucleons from weak decays 
should be excluded in the context of antideuteron pro¬ 
duction, and we thus set all particles with mean lifetimes 
above 100 fm/c to be stable. 


A. ALICE 


Deuteron and antideuteron spectra in pp minimum bias 
events at 0.9, 2.76 and 7 TeV have been measured by the 
ALICE experiment at the LHC [7]. The ALICE data 
are particularly interesting, as we here have measure¬ 
ments of both deuteron and antideuteron yields at dif¬ 
ferent energies within the same experiment. Since Monte 
Carlo antideuteron production in the coalescence model 
has shown a strong dependence on the process, there has 
been some speculation as to whether or not the coales¬ 
cence model can reproduce both the deuteron and an¬ 
tideuteron spectrum in a single experiment with the same 
value of Po- While, as can be glanced from Tables m 
and |VIIH the coalescence model yields poor fits to the 
high-energy data, the best fit values of po (and (Tq) are 
consistent between deuterons and antideuterons at the 
different energies in both Monte Carlos. The implication 
is that any future deuteron data will be very valuable for 
testing our model, or any other model of antideuteron 
formation. 

In addition to a minimum bias selection, the ALICE 
analysis uses a trigger (VOAND), which suppresses sin¬ 
gle diffractive (SD) events by requiring activity on op¬ 
posite sides of the interaction point. In order to repro¬ 
duce the results of the analysis, trigger efficiencies for 
different types of events must be taken into account. In¬ 
elastic events can be either diffractive or non-diffractive 
(ND). Since models for diffractive events, e.g. as imple¬ 
mented in Pythia 8, produce orders of magnitude fewer 
antideuterons (per event) than non-diffractive events, we 
make the approximation that only ND events will con¬ 
tribute to the antideuteron spectrum. We thus generate 
pure non-diffractive Monte Carlo events, and re-scale the 
result according to the fraction of triggered events that 
are non-diffractive, so that in this approximation predic¬ 
tions for the measured per-event spectrum can be found 


For the CLEO experiment, we set the limit at lA/c to allow 
antideuterons in T(1S) decays which is the focus of that experi¬ 
ment. 


as 


1 d^Nd 
27rAev dpTdy 


^ f 1 d-^Nd 

trig ~ 27rAfev dprdy 


(19) 


ND 


Here, A^ev is the total number of recorded/simulated 
events and Nd the number of these events with an¬ 
tideuterons. The fraction of triggered events that are 
non-diffractive is given by 


/ND.trig = 


-/VND.trig _ ENdA^ND _ SNd/nD 
A^trig Ei ’ 


( 20 ) 


where e^, W and /t, respectively, are the trigger efficiency, 
event count and fraction of the total number of events 
that are of process type i. Event counts and event frac¬ 
tions with the ‘trig’ subscript are events after trigger; the 
others are before trigger. The sum in the denominator is 
over all inelastic processes: single-, double-, central- (if 
applicable) and non-diffractive events. 

Trigger efficiencies for the different processes at 0.9, 
2.36 and 7 TeV have been estimated using Pythia 6 
and PHD JET in Ref. [H]. Using the trigger efficiencies 
and event fractions for the VOAND trigger from Tables 
5.2-5.8 in Ref. [IT], we calculate the estimated fraction 
of triggered events that are non-diffractive according to 
Eq. (20). The results are listed in Tab. jlXl In our calcu¬ 


lations of the ALICE (anti) deuteron spectra, we use the 
average value of the two Monte Carlo estimates. For the 
2.76 TeV antideuteron events, we use the event fractions 
calculated for 2.36 TeV as an estimate. 

In Fig. we show the best fits of both the coalescence 
model and our cross section based model to the ALICE 
antideuteron data at all three energies using Herwig++ 
and Pythia 8. The fits are done individually at each 
energy, and the resulting fit values can be found in Ta¬ 


bles VH and VHI For Herwig++, the slope of the spec¬ 
trum in the coalescence model is quite different from the 
slope in the experimental data, leading to very bad fits 
for the large 2.76 and 7 TeV data-sets. The cross section 
model gives a slope that is much closer to the experimen¬ 
tal result, however, the fit at 7 TeV is still quite poor for 
Herwig++. For Pythia 8, both models give better fits. 
The shape of the spectrum in the coalescence model still 
does not match the experimental data. The cross section 
based model reproduces the shape of the spectrum far 
better, and gives individual fits with y^-values that are 
consistent with the data, although visually a systematic 
overshoot at low px, and undershoot at high px seems 
to be present. This may indicate that ALICE error esti¬ 
mates are somewhat large. 

In Fig.|^we show the antiproton data from ALICE at 7 
TeV |7] compared to the spectra generated by Herwig++ 
and Pythia 8. We observe that there are small but sys¬ 
tematic differences, most significant for Herwig++. In 
terms of pT-values, these match roughly the intervals 
where the fits to the antideuteron data are poorest. Since 
the antiproton (and antineutron) spectra are the basis for 
the antideuteron spectrum, it is unreasonable to expect 
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FIG. 7: Antideuteron spectra at ALICE for the best fit values of po and (jq given in Tables [Vn| and |VIII| Top: 0.9 TeV, middle: 
2.76 TeV, bottom: 7 TeV. Left: Herwig++, right: Pythia 8. 
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Experiment 

Data points 

Best fit Po [MeV] 

xl„ 

Best fit 1/cto [barn ^ 

] Xa-Q 

d, ALICE, 0.9 TeV 

3 

228 

1.60 

5.85 

0.71 

d, ALICE, 0.9 TeV 

3 

229 

7.53 

6.15 

5.88 

d, ALICE, 2.76 TeV 

7 

199 

39.8 

4.03 

10.9 

d, ALICE, 2.76 TeV 

7 

200 

74.4 

4.00 

25.3 

d, ALICE, 7 TeV 

20 

181 

1001 

3.35 

231 

d, ALICE, 7 TeV 

20 

185 

488 

3.40 

97.4 

d, BABAR 

9 

94 

10.6 

0.63 

9.01 

d, CERN ISR 

4+4 

274 

5.15 

9.00 

5.90 

d, CLEO 

5 

130 

7.04 

0.90 

2.11 

d, LEP 

1-Pl 

152 

3.61 

1.93 

3.53 

TABLE VII: Best fit parameters for the coalescence model and the cross section model 

in Herwig++ 

Experiment 

Data points 

Best fit Po [MeV] 

Xlo 

Best fit l/(To [barn ^ 

1 Xaa 

d, ALICE, 0.9 TeV 

3 

201 

3.16 

3.58 

0.77 

d, ALICE, 0.9 TeV 

3 

201 

8.84 

3.63 

5.35 

d, ALICE, 2.76 TeV 

7 

194 

23.7 

2.93 

9.22 

d, ALICE, 2.76 TeV 

7 

196 

46.4 

2.88 

14.1 

d, ALICE, 7 TeV 

20 

194 

344 

2.63 

55.1 

d, ALICE, 7 TeV 

20 

195 

113 

2.58 

12.7 

d, BABAR 

9 

145 

16.8 

1.13 

10.1 

d, CERN ISR 

4+4 

151 

2.72 

2.08 

3.26 

d, CLEO 

5 

133 

1.16 

1.25 

1.32 

d, LEP 

H-1 

183 

3.27 

1.80 

3.55 


TABLE VIII: Best fit parameters for the coalescence model and the cross section model in Pythia 8. 


Energy //ND.trig 

Pythia 6 

PHDJET 

Average 

0.9 TeV 

0.837 

0.856 

0.847 

2.36 TeV 

0.832 

0.875 

0.854 

7 TeV 

0.831 

0.891 

0.861 


TABLE IX: Estimated fraction of minimum bias events that 
pass the ALICE VOAND trigger that are non-diffractive. 


that we can reproduce the antideuteron better than the 
progenitors. It is interesting to speculate if one could re¬ 
tune the generators so that they would better reproduce 
the antiproton data, similar to [1], and how significant 
the improvement would be for the antideuterons. How¬ 
ever, this is outside the scope of the present paper, which 
focuses on the antideuteron production model itself. 

As previously mentioned, while the coalescence model 
can give good fits to individual experiments, the best fit 
value of the coalescence momentum po varies strongly be¬ 
tween different experiments. This has lead to speculation 
that there should be some dependence on the COM en¬ 
ergy of the process in the anti deuteron production model. 
In Herwig++, (see Table VII) we see a significant energy 
dependence in the best fit values of po and l/co, with 
higher energies requiring lower values. In the coalescence 
model, the spectrum scales roughly as Pp, while in the 
cross section model it scales as I/cto- Based on this, the 



Pt [GeV] 

FIG. 8: Antiproton spectra from ALICE at 7 TeV compared 
to Monte Carlo models. 


energy dependence is significantly stronger in the coales¬ 
cence model than in the cross section model. In contrast, 
in Pythia 8, we see a much weaker energy dependence 
in both models. Here the energy dependence is some¬ 
what more significant in the cross section based model. 
While we see some sign of energy dependence in both 
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deuteron and antideuteron production, it is still unclear 
if any of this dependence can be attributed to the pro¬ 
duction models themselves, or if this is entirely an effect 
of the Monte Carlos used. 


B. BABAR 


Antideuteron production in 'T{lS,2S,iS) decays and 
non-resonant e+e“ —qq processes at y/s = 10.58 GeV 
has been measured by the BABAR Collaboration [42] . 
The latter is of particular interest for DM studies, as 
it resembles the primary annihilation process in many 
DM scenarios with a two-particle colorless (electro)weak 
initial state. It can also be directly compared to LEP 
results (see below) at -y/s = 91.2 GeV. 


The cross section based model gives better fits to the 
data than the coalescence model in both Monte Carlos, 
see Fig. and Tables VII and VIII In Herwig++, the dif¬ 
ference is rather small as the spectra from the two mod¬ 
els differ mainly at low energies, where the experimental 
data fluctuates with large errors in the two lowest energy 
bins. With Pythia 8, the cross section based model gives 
a notably better fit, as it gives a better description of the 
high energy data. 


In Pythia 8, we find the best fit parameters to the con¬ 
tinuum process to be reasonably similar to those obtained 
from T(IS') decays at CLEO (see below). This is not en¬ 
tirely unexpected, given that the two processes have very 
similar energies and a colorless initial state. However, in 
Herwig++, we find the best fit values of BABAR to lie far 
below the best fit values from CLEO. This difference in 
best fit values correspond to a factor ^ 3 in antideuteron 
production for the coalescence model, while the difference 
is around 50% in the cross section model. The explana¬ 
tion for this that is nearest at hand, is that while the 
processes have similar energies, T(IS') decays into glu- 
onic final states rather than quark final states, and that 
this brings into play differences in the two Monte Carlo 
generators. We have checked that the two Monte Carlos 
produce similar antinucleon multiplicities in each of the 
two experiments, which seems to imply that there is a 
strong process dependence in the two-(anti)baryon cor¬ 
relations from the Herwig++ cluster hadronization model 
at these energies. 

Compared to the ALICE data from proton-proton col¬ 
lisions the fitted values of pq and I/ctq are significantly 
smaller for both generators. This trend continues below 
for the other e+e^-experiments. It is difficult to deter¬ 
mine if this is somehow a process dependence that should 
be incorporated into a more complete model, or alter¬ 
natively an energy dependence, as the e+e^-experiments 
are typically at much lower energies. The latter has some 
support in the trend for better agreement at larger values 
of COM-energy seen for the LEP results. 


C. CERN ISR 


The antideuteron spectrum in pp-collisions at y/s = 
53 GeV were measured at the CERN Intersecting Storage 
Rings (ISR) at 0iab = 90° [33] and 0iab = 62.5° [33] . 

In our Monte Carlo analysis, we generate minimum 
bias events. As discussed in the ALICE analysis sec¬ 
tion, antideuterons are hardly produced in diffractive 
events, and we therefore generate purely non-diffractive 
events, and use the corresponding non-diffractive Monte 
Carlo cross sections in calculating the invariant cross sec¬ 
tion E(Pa/dp ^We note that the cross section from 
Pythia 8 is a factor ~ 22% larger than the cross sec¬ 
tion from Herwig++. Since the yield is absolute, and 
not per event, this leads to an artificial difference in the 
antideuteron yield that is not related to the event gen¬ 
eration itself, and the difference in cross section should 
therefore be kept in mind when comparing best fit pa¬ 
rameters of the two antideuteron production models. 

We perform a combined fit to the two datasets, and 
the spectra are plotted using the combined best fit values 


of Pq and ctq in Figs. 10 and II The two antideuteron 
formation models produce quite similar results in both 
Monte Carlos, and due to the large experimental errors, 
the differences in are small — the coalescence model 
giving a slightly better fit. We find the two Monte Car¬ 
los to give wildly different best fit values of po and (To: 
Herwig++ gives unusually large best fit values, whereas 
Pythia 8 gives moderately low values compared to AL¬ 
ICE. The difference in cross section between the Monte 
Carlos constitutes only a small part of this difference. 
We have checked that Pythia 8 produces a 37% higher 
multiplicity of antinucleons per event at this energy, and 
this difference is likely responsible for a sizeable part of 
the discrepancy. 

The best fit values of po and ctq differ significantly from 
the values found for the ALICE measurements in both 
Monte Carlos. In Herwig++, we see a continuation of the 
trend of increasing values with decreasing COM energies, 
and this may be another indication of an energy depen¬ 
dence stemming from the cluster hadronization model. 
While we saw a similar tendency in the Pythia 8 ALICE 
results, the ISR results do not support the hypothesis of a 
systematic COM energy dependence in this Monte Carlo. 


D. CLEO 


Antideuteron production 
measured at CLEO [35] . 
the CLEO data are shown 
Pythia 8 have similar best 


in T(IS') decays has been 
The best fit spectra for 
in Fig. 12 Herwig++ and 


fit values of po in the coa- 


In our previous work [4], we used the larger experimentally mea- 
sured total inelastic cross section, leading to an overestimation 
of the antideuteron yield. 
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FIG. 9: Antideuteron spectra at BABAR for the best fit values of po and (Tq given in Tables m and |VIII[ Left: Herwig++, 
right: Pythia 8. 
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FIG. 10: Antideuteron spectra at ISR, generated using Herwig++ with the best fit values of po and cro given in Tab. IVlI] Left: 
^lab = 90°, right: Slab = 62.5°. 


lescence model, while the best fit (Tq differs quite signif¬ 
icantly in the cross section based model. In Herwig++, 
the cross section based model reproduces the shape of the 
spectrum significantly better than the coalescence model, 
and thus gives a better fit. In Pythia 8, the coalescence 
model gives a better fit to the low energy data, while the 
cross section based model gives a better fit to the high 
energy data. As a result the two models give very similar 
fits; the coalescence model having a slightly lower 


E. LEP 

Antideuteron searches in e“'"e“-collisions were per¬ 
formed by the ALEPH [1^ and OPAL [17] experiments at 
LEP. Both collaborations studied antideuteron multiplic¬ 
ities in hadronic events at the Z-resonance. ALEPH ob¬ 
served (5.9± 1.8±0.5) X 10“® antideuterons per hadronic 
event in the momentum range 0.62 < p < 1.03 GeV and 
angular range |cos0| < 0.95, where the errors are sta¬ 
tistical and systematical, respectively. In OPAL, how¬ 
ever, no antideuteron candidates were observed in the 
antideuteron momentum range 0.35 < p < 1.1 GeV. In 
previous works, only the ALEPH result has been used 


for calibration of the the coalescence momentum, but the 
negative OPAL result should also be taken into account. 
As the expected number of signal events in the two exper¬ 
iments are of the same magnitude, the non-observation 
of antideuterons in OPAL might be an indication that 
the ALEPH result suffers from an upwards fluctuation. 
Performing a combined fit of the two experiments will 
yield a lower best fit coalescence momentum than previ¬ 
ous fits based on the ALEPH data alone. 

In order to calculate the for OPAL, we first estimate 
the expected number of signal events s at OPAL by 


S — 


( 21 ) 


where e = 0.234 is the given detection efficiency, N^y = 
1.64 X 10® is the number of events in the OPAL analysis, 
and nj is the Monte Carlo prediction for the number 
of antideuterons per event, e and njjyjQ are in reality 
energy dependent quantities, but only the average value 
of the detection efficiency is available. Using the fact that 
no antideuteron candidates were observed by OPAL, and 
assuming Poissonian uncertainty a = y/s for the expected 
number of events, the x^ is then given by 


2 (fVobs - s)^ 

XOPAL — o -“ “ '®- 


(22) 
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FIG. 11: Antideuteron spectra at ISR, generated using Pythia 8 with the best £t values of po and ao given in Tab. |VIII| Left: 
^lab = 90°, right: Slab = 62.5°. 



FIG. 12: Antideuteron spectra at CLEO for the best fit values of po and ao given in Tab. ivnl Left: Herwig++, right: Pythia 8. 


F. Combined fits 

For the purpose of calculating the cosmic ray an¬ 
tideuteron flux from dark matter, the free parameters 
Po and cTo have to be calibrated based on fits to the 
previously discussed experimental data. This calibra¬ 
tion should be done separately for each Monte Carlo, 
as the best fit values generally differ between Monte 
Carlos, which should be clear from the above, e.g. due 
to differences in primary antinucleon spectra and the 
(anti)nucleon correlations. It is also clear that no pa¬ 
rameter values exist that give good simultaneous fits to 
all experiments, and this can be seen quantitatively in 
Table ra combined fits to all experiments yield y^/d.o.f 
ranging from 10 to a whopping 54. 

In order to get sensible results, it is necessary to re¬ 
strict the fits to reasonably self-consistent subsets of ex¬ 
periments. However, it is not a priori clear which exper¬ 
iments should be included in the fits. In previous work, 
the choice of po for the coalescence model has often been 
based on a fit to the ALEPH data alone, as LEP events 
are similar to DM annihilation events in typical DM mod¬ 
els. However, as discussed earlier, the OPAL experiment 
at LEP did not observe any antideuterons in a similar 
range of energies. In fact, even the ALEPH data alone 


does not constitute more than a 3cr observation of an¬ 
tideuterons. The problem of relying on a single data 
point for calibration has also been discussed extensively 
in the past, e.g. see [9j. 

We will here divide the experiments into two groups: 
experiments with colored initial states (ALICE, ISR), 
and experiments with colorless initial states (BABAR, 
CLEO, LEP), and consider separate fits to these two sets. 
While the antideuteron formation process is in both mod¬ 
els assumed to be agnostic to the nature of the hard pro¬ 
cess, it is not unlikely that the outcomes of the hadroniza- 
tion models of the Monte Carlos are sensitive to differ¬ 
ences in the underlying physics between these two classes 
of processes, and thus originate a difference in the fitted 
value. As dark matter annihilations have colorless ini¬ 
tial states, the colorless set is likely the most relevant for 
calculating the cosmic ray antideuteron flux from dark 
matter. 

Best fit values for the two sets of experiments can be 
seen in Table In all cases, the cross section based 
model gives a considerably better combined fit than the 
coalescence model. In Herwig++, the fits are still rather 
bad for both datasets with either model. This is not en¬ 
tirely unexpected; in the colorless set, the BABAR data 
prefers much lower values of the free parameters than the 
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Monte Carlo Experiments 

Data points 

Best fit Po [MeV] 

Xpn 

Best fit 1/cro [barn 


Herwig++ ALICE (d), ISR 

38 

187 

646 

3.50 

196 

BABAR, CLEO, LEP 

16 

96 

73.6 

0.68 

29.2 

All experiments 

54 

123 

2859 

1.43 

2146 

Pythia 8 ALICE (d), ISR 

38 

193 

255 

2.63 

58.2 

BABAR, CLEO, LEP 

16 

140 

30.5 

1.18 

16.7 

All experiments 

54 

174 

888 

2.13 

510 


TABLE X: Combined best fit parameters for the coalescence model and the cross section model for different selections of 
experimental data. 


other experiments, thus giving a bad simultaneous fit. In 
the set with colored initial states, the poor individual fits 
of Herwig++ to the ALICE data alone are enough to give 
a bad combined fit, and the large spread in the individual 
best fit parameters further worsens the result. 

In Pythia 8, the coalescence model gives relatively 
poor fits to both sets. The cross section based model, 
on the other hand, gives a good fit to the set with color¬ 
less initial states with = 16.7 for 15 degrees of freedom 
(d.o.f), and gives a decent fit with = 58.2 for 37 d.o.f 
to the set with colored initial states. 


V. DARK MATTER SPECTRA 

We will here compare the antideuteron spectra at 
Earth coming from a generic dark matter candidate anni¬ 
hilating into bb and W^W~ in the coalescence and cross 
section based models, in order to see the impact of the 
new model on the, in principle, measurable spectrum. We 
will be using the best fit values of po and ctq discussed in 
Sec. |IVF[ and we will consider dark matter candidates 
with masses of 100, 500 and 1000 GeV. 


A. Propagation of antideuterons 

Antideuterons, being charged particles, do not prop¬ 
agate through our galaxy in straight lines, but are de¬ 
flected in the turbulent Galactic magnetic fields. This 
leads to a random walk behaviour, which can be well 
described using a diffusion model. The most commonly 
used model is the so-called two-zone diffusion model - a 
cylindrical model consisting of a magnetic halo region of 
radius R = 20 kpc and half-height L, where charged par¬ 
ticles diffuse freely; and a thin gaseous disk of the same 
radius and a half-height of h = 100 pc, where scattering 
and annihilation on interstellar matter can additionally 
take place. While R and h are set by the size of the 
observed Galactic disk, the half-height of the magnetic 
halo, L, is a free parameter. 

For antideuterons, energy redistribution terms and 
non-annihilating inelastic scattering only constitute mi¬ 
nor corrections, and are typically neglected, as they will 
be here. Under the assumption of steady state condi¬ 
tions, the diffusion equation describing this model is then 


given by 

- D{T)W^f + ^(sign(z)/W) = Q - 2W(z)r,„,(T)/ , 

(23) 

where f{x,T) = dN^/dT is the number density of an¬ 
tideuterons per unit kinetic energy T, D{T) = Dq[3TZ^ 
is the (spatial) diffusion coefficient, Q is the source term 
from dark matter annihilations, W is the velocity of a 
convective wind perpendicular to the Galactic disk, z is 
the vertical coordinate, (3 = v/c is the antideuteron ve¬ 
locity, and TZ is the antideuteron rigidity in units of GV. 
(5, Dq, and 14 are here free parameters of the model. 

The annihilation rate, Tann, of antideuterons on inter¬ 
stellar gas in the Galactic disk is given by 

ran„(T) = {riH + 4inHe)t'a|“(r), (24) 

where uh ~ 1 cm“^ and « 0.07n/f are the respective 
number densities of hydrogen and helium nuclei in the 
disk. The factor 43 here accounts for the difference in 
annihilation cross section between H and He, under the 
assumption of simple geometrical scaling. We estimate 
the annihilation cross section using 

- n™'’"°"-“"(r), (25) 

where = a{dp —)■ dX) is the component of 

the inelastic cross section where the antideuteron sur¬ 
vives the interaction. Data on these cross sections are 
sparse, and it is therefore necessary to make approxima¬ 
tions based on re-scaling of pp data and use of charge 
conjugate processes. This has recently been discussed in 
detail in Ref. [33], and we will here adopt their fits to 
experimental data. 

The source term Q is for the case of annihilating dark 
matter given by 


Q{r,T) 


lp2(f) 


2 


^(cru). 


dm 

_ d 

dT ’ 


(26) 


where p{r) is the dark matter density, mpM is its mass, 
and {av)i is the thermally averaged dark matter annihila¬ 
tion cross section for channel i. For the dark matter halo 
profile, we chose the Navarro-Frenk-White (NFW) [48] 
profile. 


p(r) 


Po 

(r/rs) [1 + ir/rs)f' 


(27) 
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Model 

L in kpc 

S 

Do in kpc'^ Myr ^ 

Vc in kms ^ 

max 

15 

0.46 

0.0765 

5 

med 

4 

0.7 

0.0112 

12 

min 

1 

0.85 

0.0016 

13.5 


TABLE XI: Propagation parameters for the max, med and 
min models. 


which has been shown to be in good agreement with 
the results of N-body halo formation simulations. For 
the free parameters in the NFW-profile we use po = 
0.26 GeV/cm^ and rs = 20 kpc. 

For the free parameters of the diffusion model, it has 
been common in the literature to use the three sets of 
values given in Ref. |40j . that yield maximal, median 
and minimal antiproton fluxes from dark matter anni¬ 
hilation, while remaining compatible with the observed 
B/C ratio. These parameter sets are labeled ‘max’, ‘med’ 
and ‘min’ respectively, and their values are listed in Ta¬ 
ble IXIl The ‘max’ and ‘min’ models are often used to 
estimate the uncertainty band from propagation, but as 
these are are physically extreme models, the resulting un¬ 
certainty band is likely overly conservative. Indeed, the 
‘min’ model has recently been excluded by cosmic ray 
positron data m- Propagation uncertainty has been 
thoroughly discussed in the literature, and is not the fo¬ 
cus of this article. We therefore restrict our propagation 
calculation to the ‘med’ model. 

The diffusion equation (23) can be solved semi- 
analytically m, and for annihilating dark matter, the 
expression for antideuteron flux near Earth is 




Po \ 
muM/ 


R{T)- 


dN^ 
2 dT 


(28) 


where 


OO 

RiT) = Y,Jo[cJf) 


exp - 


VaL 


VniL) 


yn{z) = 


jncn)R^ Jo 

V,{Z-z) 


2 K J A„sinh(S'„L/2)’ 
(29) 




exp 


2D 


sinh 


Sn{Z - z)\ f p{r,z) 


) (^) }■ 

(30) 


= 2/ir,„„ + P, + DSr. coth(^„L/2), (31) 

and 

The particle physics of dark matter annihilation and 
the astrophysics of the propagation are here neatly sep¬ 
arated — the astrophysics of the propagation is con¬ 
tained within the propagation function i?(T), which is 


completely independent of the particle physics of the an¬ 
nihilation process. This function can thus, independently 
of the dark matter model in question, be tabulated for 
a given halo and set of diffusion model parameters, and 
later applied to the propagation of antideuterons from 
any model of (symmetric) dark matter annihilation. 

Solar modulations is taken into account using the force 
field approximation |52j . shifting the kinetic energy of 
the particles from T to a kinetic energy near Earth of 
7® = r — |Ze|(()Fisk, where the so-called Fisk potential 
</'Fisk = 0.5 GV is an effective potential that parametrizes 
the energy loss from the solar wind. The corresponding 
antideuteron flux near Earth is then finally given by 




2 m jT -h r2 


(33) 


More realistic modeling, as well as an estimation of un¬ 
certainties due to solar modulation has been discussed in 
detail in Ref. [6]. 


B. Antideuteron flux near Earth 


We generate events for dark matter annihilations into 
bb and W~^W~ final states for DM masses of 100 GeV, 
500 GeV and I TeV using Herwig++ and Pythia 8. For 
the antideuteron formation we use the two sets of best 
fit values of po and g p, based on pp and e+e^-data, as 
discussed in Sec. IV F We assume 100% branching ratios 
into the given channels, and use the canonical value of 
{av) = 3 X 10“^® cm^s“^ for the thermally averaged dark 
matter annihilation cross section. 

In Figs. [T3|and[T^we, respectively, show the expected 
antideuteron fluxes after propagation from Herwig++ and 
Pythia 8. The figures show the fluxes as a function of 
the kinetic energy per nucleon of the antideuteron in both 
the coalescence model and the cross section based model, 
using the calibration of uo and po against experiments 
with colorless initial states (BABAR, GLEO, LEP). The 
bands indicate the statistical uncertainty in our event 
generation, and the shaded regions at the top indicate 
the most recent values for the expected sensitivities of 
the AMS-02 and GAPS experiments [53]. 

When comparing the predicted fluxes from the two 
models, one should keep in mind that the relative normal¬ 
ization of the fluxes is not fixed, but determined by the 
calibration of the free parameters po and ctq. The shapes 
of the spectra are, however, more or less independent of 
the calibration, and comparing the shapes thus gives a 
more reliable picture of the difference between the mod¬ 
els. In particular, the differences between the two models 
appear larger in Herwig++ than in Pythia 8, but this is 
largely an effect of the calibration. Figure shows the 
Herwig++ result using the calibration against colored ini¬ 
tial states (ALICE, ISR), and we see that the difference 
between the two models is considerably smaller here due 
to less of a difference in normalization. In Pythia 8, the 
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T [GeV/n] T [GeV/n] 


FIG. 13: Antideuteron spectra at Earth from dark matter annihilation into bb (left) and W'^W~ (right), calculated using 
Herwig++ with po = 96 MeV in the coalescence model, and I/ctq = 0.68 barn“^ in the cross section based model. The dashed 
line shows the expected astrophysical background calculated in Ref. | 54| . 




T [GeV/n] T [GeV/n] 


FIG. 14: Antideuteron spectra at Earth from dark matter annihilation into bb (left) and W'^W~ (right), calculated using 
Pythia 8 with po = 140 MeV in the coalescence model, and 1/cro = 1.18 barn“^ in the cross section based model. The dashed 
line shows the expected astrophysical background calculated in Ref. | 54| . 


difference in normalization between the models is simi¬ 
lar in the colored and colorless calibrations. We therefore 
leave out the plot for the colored calibration in Pythia 8. 

The differences in the shapes of the spectra between 
the two models appear to be similar in the two Monte 
Carlos. In the bb channel, we see a consistent qualitative 
difference between the models across all DM masses: the 
cross section based model predicts a softer antideuteron 
spectrum, with a more rapid falloff at high energies. The 
same can be seen in the 100 GeV dark matter W'^W~ 
final-state. This leads to an enhanced flux in the low 
energy range relevant for AMS-02 and GAPS, where the 
background is expected to be small. With the values of 
Po and CTo used here, the predictions for the flux from 
the two models typically differ by a factor 1.5-2 in the 
experimentally relevant energy ranges. 


For the higher masses, the situation is less clear for 
the W~^W~ final-state due to the statistical uncertainty 
from limited statistics in the Monte Carlo event genera¬ 
tion. The two models seem to predict similar slopes at 
low energies, but the cross section model shows signs of 
a steeper falloff at high energies. We see that the cross 
section model consistently predicts a higher flux at the 
peak than the coalescence model. This leads to a possi¬ 
bly enhanced flux compared to the coalescence model in 
the multi-GeV kinetic energy region where the AMS-02 
experiment has some sensitivity. 
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FIG. 15: Antideuteron spectra at Earth from dark matter annihilation into bb (left) and W'^W~ (right), calculated using 
Herwig++ with po = 187 MeV in the coalescence model, and l/no = 3.50 barn“^ in the cross section based model. The dashed 
line shows the expected astrophysical background calculated in Ref. | 54| . 


VI. CONCLUSIONS 

We have proposed a new model for describing the 
formation of antideuterons in high energy events. Our 
model is based on the experimentally measured cross sec¬ 
tions for nucleon capture processes, and is — in contrast 
to the state-of-the-art coalescence model — capable of 
describing recent deuteron and antideuteron data from 
the ALICE experiment at the LHC. 

The physical interpretation of the antideuteron forma¬ 
tion process differs significantly between our model and 
the coalescence model. In the coalescence model, an¬ 
tideuteron formation is described by slow nucleon cap¬ 
ture, whereas in our model, antideuterons are primar¬ 
ily produced through resonant processes with the delta- 
resonance, which peaks for COM momentum differences 
near 1 GeV. Moreover, while the coalescence model 
strictly describes a pn capture process, our model pre¬ 
dicts similar antideuteron contributions from pp and fin 
processes. 

We have compared the predictions of our model to the 
coalescence model for several different experiments, and 
find our model to give comparable or better descriptions 
of the data in all experiments; the difference being most 
significant for the ALICE experiment, where the coales¬ 
cence model fails to give a satisfactory description. For 
the purpose of dark matter indirect detection, we perform 
fits of the free parameters of the models against two sets 
of experimental data, divided into experiments with col¬ 
ored and colorless initial states. We find our model to 


give consistently better simultaneous fits to the exper¬ 
imental data in both Herwig++ and Pythia 8, and in 
Pythia 8, the fits for our cross section based model give 
X^-values that indicate the model can describe the data 
successfully. 

Comparing the predicted antideuteron fluxes from 
dark matter annihilation in the two models, with a se¬ 
lection of different dark matter masses and different final 
states, we hnd that our model produces softer spectra 
than the coalescence model, thus giving an enhanced an¬ 
tideuteron flux in the low kinetic energy range relevant 
for current and planned experiments. 
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